#Created by Lu Shen, Oct 20
linear.detrend=function(temp){
	xx=1:length(temp)
	MAT=na.omit(cbind(xx,temp))
	spdata=array(NA,c(length(xx)))
	if (length(MAT)>=8){
		fit<-lm(MAT[,2]~MAT[,1])
		yy=xx*fit$coefficients[2]+fit$coefficients[1]
		spdata=temp-yy
	}
	return(spdata)
}

mov.detrend=function(temp,len=7){
	  ap=mov.avg.wgt(temp, rep(1,len))	  
	  return(as.numeric(temp-ap))
}

find.cor2=function(spdata,index,p_value){
	new=array(NA,dim(spdata))
	for (i in 1:dim(new)[1]){
		for (j in 1:dim(new)[2]){
			if (sum(!is.na(spdata[i,j,]))>=5){
			new[i,j,]=index
			}
		}
	}
	corr=find.cor(spdata,new,plim=p_value)
	return(corr)
}

find.cor3=function(spdata,y){
	correlation=array(NA,c(dim(spdata)[1],dim(spdata)[2]))
	p_values=array(NA,c(dim(spdata)[1],dim(spdata)[2]))	
	for(ii in 1:dim(spdata)[1]){
		for(jj in 1:dim(spdata)[2]){
			if(sum(!is.na(spdata[ii,jj,]))>=5){
			fit<-cor.test(spdata[ii,jj,],y)
			correlation[ii,jj]=fit$estimate
			p_values[ii,jj]=fit$p.value
			}
		}
	}
	return(list(corr=correlation,p=p_values))
}

test_siginificance=function(spdata,mu,plim){
significance=array(0,c(dim(spdata)[1],dim(spdata)[2]))
for (i in 1:dim(spdata)[1]){
	for (j in 1:dim(spdata)[2]){
		y=spdata[i,j,]
		if (sum(!is.na(y))>5){
			ap=t.test(y,mu=mu, alternative="two.sided")$p.value
			if (ap<plim) significance[i,j]=1
		}
	}	
}
return(significance)	
}

read_US=function(name,pressure,start_month,end_month,xlim=c(210,290),ylim=c(10,70)){
date=array(NA,c(792,2))
for (i in 1:792){
  if (i%%12>0) {date[i,1]=1948+floor(i/12);date[i,2]=i%%12}
  if (i%%12==0){date[i,1]=1948+floor(i/12)-1;date[i,2]=12}
}	
datadir='./Data/pressure/'
datafile = nc_open(paste(datadir,name,sep=''))
lon = ncvar_get(datafile, varid='lon')
lat = ncvar_get(datafile, varid='lat')
ind1=which(lon>=xlim[1] & lon<=xlim[2])
ind2=which(lat>=ylim[1] & lat<=ylim[2])
ind4=which(date[,1]>=1948 & date[,1]<=2013)
lon=lon[ind1]
lat=lat[ind2]
date=date[ind4,]
level = ncvar_get(datafile,varid='level')
ind3=which(level== pressure)
origin_data = ncvar_get(datafile,start=c(ind1[1],ind2[1],ind3,ind4[1]), count=c(length(ind1),length(ind2),1,length(ind4)))
origin =origin_data[,length(lat):1,]
lat=rev(lat)
nc_close(datafile)
temp=array(NA,c(length(lon),length(lat),66))
for (year in 1948:2013){
  if(start_month<=end_month){
	  ind=(date[,1]==year & date[,2]>=start_month & date[,2]<=end_month)
  } else {
  	 ind=((date[,1]==(year-1) & date[,2]>=start_month) | (date[,1]==year & date[,2]<=end_month))
  }
  
  if(sum(ind)>1) temp[,,year-1947]=apply(origin[,,ind],c(1,2),mean,na.rm=TRUE)
  if(sum(ind)==1) temp[,,year-1947]=origin[,,ind]
  
}
return(list('data'=temp,'lon'=lon,'lat'=lat))
}


read_surface=function(name,start_month,end_month,xlim=c(210,290),ylim=c(10,70)){
date=array(NA,c(792,2))
for (i in 1:792){
  if (i%%12>0) {date[i,1]=1948+floor(i/12);date[i,2]=i%%12}
  if (i%%12==0){date[i,1]=1948+floor(i/12)-1;date[i,2]=12}
}
datadir='./Data/surface/'	
datafile = nc_open(paste(datadir,name,sep=''))
lon = ncvar_get(datafile, varid='lon')
lat = ncvar_get(datafile, varid='lat')
ind1=which(lon>=xlim[1] & lon<=xlim[2])
ind2=which(lat>=ylim[1] & lat<=ylim[2])
ind4=which(date[,1]>=1948 & date[,1]<=2013)
lon=lon[ind1]
lat=lat[ind2]
date=date[ind4,]
origin_data = ncvar_get(datafile,start=c(ind1[1],ind2[1],ind4[1]), count=c(length(ind1),length(ind2),length(ind4)))
origin =origin_data[,length(lat):1,]
lat=rev(lat)
nc_close(datafile)

temp=array(NA,c(length(lon),length(lat),66))
for (year in 1948:2013){
  if(start_month<=end_month){
	  ind=(date[,1]==year & date[,2]>=start_month & date[,2]<=end_month)
  } else {
  	 ind=((date[,1]==(year-1) & date[,2]>=start_month) | (date[,1]==year & date[,2]<=end_month))
  }

  if(sum(ind)>1) temp[,,year-1947]=apply(origin[,,ind],c(1,2),mean,na.rm=TRUE)
  if(sum(ind)==1) temp[,,year-1947]=origin[,,ind]

  # temp[,,year-1947]=apply(origin[,,ind],c(1,2),mean,na.rm=TRUE)
}
return(list('data'=temp,'lon'=lon,'lat'=lat))
}

read_SST=function(name,start_month,end_month,xlim=c(210,330),ylim=c(10,70)){
date=array(NA,c(1933,2))
for (i in 1:1933){
  if (i%%12>0) {date[i,1]=1854+floor(i/12);date[i,2]=i%%12}
  if (i%%12==0){date[i,1]=1854+floor(i/12)-1;date[i,2]=12}
}
datadir='./Data/SST/'
datafile = nc_open(paste(datadir,name,sep=''))
lon = ncvar_get(datafile, varid='lon')
lat = ncvar_get(datafile, varid='lat')
time.vec = ncvar_get(datafile, varid='time')
ind1=which(lon>=xlim[1] & lon<=xlim[2])
ind2=which(lat>=ylim[1] & lat<=ylim[2])
lon=lon[ind1]
lat=lat[ind2]
ind3=which(date[,1]>=1948 & date[,1]<=2013)
origin_data = ncvar_get(datafile,start=c(ind1[1],ind2[1], ind3[1]),count=c(length(ind1),length(ind2),length(ind3)))
nc_close(datafile)
new_data=origin_data[,length(lat):1,]
temp=new_data
lat=rev(lat)
date=date[ind3,]
SST=new_data

year_SST=array(NA,c(dim(SST)[1],dim(SST)[2],66))
for (year in 1948:2013){
  if(start_month<=end_month){
	  ind=(date[,1]==year & date[,2]>=start_month & date[,2]<=end_month)
  } else {
  	 ind=((date[,1]==(year-1) & date[,2]>=start_month) | (date[,1]==year & date[,2]<=end_month))
  }
  if(sum(ind)>1) year_SST[,,year-1947]=apply(SST[,,ind],c(1,2),mean,na.rm=TRUE)
  if(sum(ind)==1) year_SST[,,year-1947]= SST[,,ind]  
	# year_SST[,,year-1947]=apply(SST[,,ind],c(1,2),mean,na.rm=TRUE)	
}
return(list('data'= year_SST,'lon'=lon,'lat'=lat))
}

read_20c=function(name,start_month,end_month,xlim=c(220,320),ylim=c(10,70)){
time=seq(as.Date('1871-01-01'),as.Date('2012-12-31'),by='month')
date=cbind(as.numeric(format(time,"%Y")),as.numeric(format(time,"%m")))	
datadir='./Data/20Century/'
#--read the data---------
datafile = nc_open(paste(datadir,name,sep=''))
lon = ncvar_get(datafile, varid='lon')
lat = ncvar_get(datafile, varid='lat')
ind1=which(lon>=xlim[1] & lon<=xlim[2])
ind2=which(lat>=ylim[1] & lat<=ylim[2])
ind4=which(date[,1]>=1871 & date[,1]<=2012)
lon=lon[ind1]
lat=lat[ind2]
date=date[ind4,]
if(datafile$ndims==5){
origin_data = ncvar_get(datafile,start=c(ind1[1],ind2[1],1,ind4[1]), count=c(length(ind1),length(ind2),1,length(ind4)))
} else {
origin_data = ncvar_get(datafile,start=c(ind1[1],ind2[1],ind4[1]), count=c(length(ind1),length(ind2),length(ind4)))	
}
origin =origin_data[,length(lat):1,]
lat=rev(lat)
nc_close(datafile)

temp=array(NA,c(length(lon),length(lat),2012-1871+1))
for (year in 1871:2012){
  if(start_month<=end_month){
	  ind=(date[,1]==year & date[,2]>=start_month & date[,2]<=end_month)
  } else {
  	 ind=((date[,1]==(year-1) & date[,2]>=start_month) | (date[,1]==year & date[,2]<=end_month))
  }
  temp[,,year-1870]=apply(origin[,,ind],c(1,2),mean,na.rm=TRUE)
}

return(list('data'=temp,'lon'=lon,'lat'=lat,'date'=1871:2012))
}

test_significance=function(spdata,cr, lon.frame=lon.frame,lat.frame=lat.frame){
	ind=(abs(spdata)> cr)
	result=matrix(0,0,2)
	for (i in 1:length(lon.frame)){
		for (j in 1:length(lat.frame)){
			if (!is.na(ind[i,j])){
				if (ind[i,j]){
					temp=c(lon.frame[i],lat.frame[j])
					result=rbind(result,temp)
				}
			}
		}
	}
	return(result)
}


cal.area=function(lon,lat){
SST_area=array(NA,c(length(lon),length(lat)))
delta1=diff(lon)[1]/2;delta2=diff(lat)[1]/2
lon0=lon[1];
if(lon0>=180) lon0=lon0-360
	for (j in 1:length(lat)){
	 x1=lon0;x2=lat[j]
     p <- rbind(c(x1+delta1,x2+delta2),c(x1+delta1,x2-delta2),c(x1-delta1,x2-delta2),c(x1-delta1,x2+ delta2)) 	
     SST_area[,j]=areaPolygon(p)
	}
	return(SST_area)
}

cal.area_mean=function(spdata,area,ind1,ind2){
	y=array(NA,dim(spdata)[3])
	for(k in 1:length(y)){
		y[k]=sum(spdata[ind1,ind2,k]* area[ind1,ind2],na.rm=TRUE)/sum(area[ind1,ind2],na.rm=TRUE)
	}
	return(y)
}
